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Abstract 

This  paper  is  concerned  with  developing  methods  for  the  propagation  of  linear  waves 
in  several  space  dimensions.  The  methods  are  time-reversible,  and  hence  free  from  nu¬ 
merical  dissipation.  They  are  based  on  bicharacteristic  forms  of  the  governing  equations, 
and  are  made  possible  by  adopting  forms  of  staggered  storage  that  depend  on  the  pre¬ 
cise  equations  under  consideration.  Analysis  is  presented  for  the  equations  of  acoustics, 
electromagnetics,  and  elastodynamics. 
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1  Introduction 


In  Computational  Physics  we  solve  numericaDy  some  discrete  replacement  of  the  governing 
equations,  often  because  those  equations  are  nonlinear,  and  hence  not  amenable  to  analytical 
treatment.  There  is  a  comparative  neglect  of  linear  problems,  except  as  preliminary  models 
on  which  to  illustrate  algorithmic  principles.  There  are  however  many  linear  problems 
that  are  technically  important  in  their  own  right,  and  which  demand  huge  computational 
resources,  due  to  the  neccessity  of  representing  very  complicated  solutions.  In  this  paper  we 
shall  be  concerned  with  problems  in  which  low  amplitude  wave  propagation  takes  place  over 
distances  typified  by  large  multiples  of  the  wavelength.  A  few  of  the  areas  and  applications 
where  such  problems  arise  are 

1.  Acoustics  and  Aeroacoustics,  for  noise  abatement 

2.  Electromagnetics,  for  microcircuit  design, 

3.  Elastodynamics,  for  nondestructive  testing, 

4.  Seismology,  for  oil  exploration, 

5.  Medical  Imaging,  for  accurate  diagnosis, 

6.  Hyperthermia,  for  noninvasive  surgery. 

In  some  of  these  case?,  the  wavelengthss  are  required  to  be  small  compared  to  the 
obstacles  they  travel  round,  to  enable  accurate  measurements,  but  if  the  wavelengths  are 
not  so  small  that  geometrical  optics  can  be  applied,  we  have  an  ‘intermediate  regime’  that  is 
computationally  the  most  demanding.  In  other  cases,  this  regime  arises  by  accident  rather 
than  design.  Either  way,  there  is  a  need  to  propagate  smaU-amplitude  waves  over  tens  or 
hundreds  of  wavelengths  with  minimal  dispersion  and  damping. 

The  view  taken  here  will  be  that,  although  both  kinds  of  error  are  important,  it  is  often 
damping  that  is  the  more  damaging  to  practical  calculations.  A  wave  that  is  excessively 
damped  simply  disappears,  and  there  is  no  solution  to  see.  Dispersion  leaves  the  wave 
intact,  but  in  the  wrong  place.  If  all  waves  were  to  propagate  at  the  same  wrong  speed,  we 
would  have  a  good  answer  to  a  question  that  differs  from  the  one  intended;  our  waves  would 
propagate  in  the  wrong  medium.  However,  since  the  data  we  have  concerning  that  medium 
is  usually  inexact,  this  may  not  be  of  great  concern.  In  this  paper,  we  consider  schemes  that 
by  construction  are  completely  free  from  dissipation  and  which  maintain  dispersion  errors 
within  some  acceptable  limits. 

The  schemes  studied  here  are  leapfrog  methods  (which  are  inherently  free  of  dissipation) 
but  combined  with  upwind  bias.  Such  schemes  were  first  described  by  Iserles[l]  for  one¬ 
dimensional  advection.  Compared  with  classical  leapfrog  methods  they  have  more  compact 
stencils,  leading  to  improved  accuracy  and  more  natural  handling  of  boundaries.  An  ac¬ 
count  of  the  one-dimensional  methods  is  provided  in  Section  2,  which  provides  a  summary, 
commentary,  and  extension  of  [1]. 

The  present  paper  gives  the  first  use  of  these  methods  in  higher  dimensions,  and  the 
generalisation  involves  some  new  ideas.  For  any  characteristic-based  scheme  there  is  a  need 
to  reconcile  computational  convenience  and  physical  motivation.  These  issues  are  discussed 
in  Section  3.  In  Section  4  it  is  shown  that  one  solution  is  to  introduce  staggered  storage, 
in  a  manner  dictated  by  the  coupling  of  unknowns  in  the  given  equation  set.  Particular 
examples  discussed  are  the  equations  of  acoustics,  elastodynamics,  and  electromagnetics. 
Computational  results  are  presented  in  Section  5. 
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Figure  1:  Four  symmetric  stencils  for  wave  propagation.  '''  J-'.ssical  leapfrog,  (b)  Upwind 
leapfrog,  (c)  Upwind  leapfrog  extended  in  space,  (d)  Up.  .e<  "og  extended  in  time. 


2  Upwind  Leapfrog  Schemes  in  One  Dimension 

2.1  Symmetry  and  Dissipation 

To  solve  any  wave  propagation  problem  without  introducing  dissipation,  it  is  obviously 
necessary  that  the  scheme  employed  be  reversible  in  time.  This  implies  that  the  stencD  on 
which  the  scheme  is  supported  has  central  symmetry,  as  is  the  case  for  the  four  stencils 
shown  in  Figure  1.  Example  (a)  is  the  stencil  for  the  classical  leapfrog  scheme.  This  has 
symmetry  both  in  space  (left-right)  and  in  time  (forward-backward).  However,  the  weaker 
(point)  symmetry  of  the  other  examples  is  also  adequate  to  ensure  reversibility.  A  scheme 
without  dissipation  can  be  used  to  integrate  from  given  data  at  t  =  0  to  a  solution  at  t  =  T, 
and  then,  with  time  reversed,  integrated  back  to  t  =  0  to  recover  the  data  exactly  (apart 
from  roundoff  errors).  During  the  reverse  integration  the  direction  of  wave  propagation 
is  also  reversed;  hence  the  requirement  is  for  the  stencil  to  remain  invariant  under  a  180° 
rotation. 

It  is  also  neccessary  for  reversibility  that  the  coefficients  of  the  scheme  be  symmetrical. 
If  the  scheme,  for  some  scalar  problem  in  which  the  unknown  is  u  is  written 

^CpWp  =  0 


then  the  condition 

Cpi  +  Cp2  =  0 

where  pl,p2  are  any  pair  of  points  that  exchange  places  under  the  rotation,  is  neccessary 
and  sufficient  to  ensure  reversibility. 

The  modified  equation  of  such  a  scheme  contains  no  even  order  derivatives;  as  an  ap¬ 
proximation  to  a  first-order  PDE  it  therefore  has  even  order  of  accuracy. 
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2.2  Errors  and  Stability  Bounds 
2.2.1  Stencil  la 

On  the  stencil  of  Figure  1(a)  the  only  consistent  schenae  for  linear  advection 

uj  +  aux  =  0 


is  the  regular  leapfrog  method 

=  (1) 

where  the  Courant  number  u  =  aAt/Ax. 

For  this  and  other  reversible  methods  the  dispersion  relationship  can  be  found  by  as¬ 
suming  solutions  of  the  form 

uj  =  exp  t(n<f>  -  jd)  (2) 

Given  a  real  value  for  the  Fourier  angle  d,  finding  a  real  value  for  the  phase  change  <t> 
will  indicate  that  the  method  is  dissipationless  as  assumed.  A  complex  value  of  <f>  means 
that  exp  icf)  is  not  on  the  unit  circle,  and  hence  that  there  are  amplification  factors  ai ,  02 
exhibiting  growth  or  decay.  Time  reversibility  implies  that  0102  =  1,  and  so  the  scheme  has 
an  unstable  mode  in  either  direction. 

Substituting  (2)  into  (1)  gives  the  dispersion  relationship  as 

sin  =  t'  sin  0 


and  it  is  evident  that  the  scheme  is  unstable  for  some  0  if  1/  >  1. 

For  1/  <  1  the  only  error  is  in  phase  speed.  A  signal  such  as  (2)  travels  with  speed 


(f>  Ax 


We  therefore  define  the  phase  speed  error  as 


£r,  = 


Ax 

At 


Ar 

At 


(3) 

(4) 


Here  the  ‘grid  speed’  Ax /At  has  been  used  to  make  the  error  non-dimensional.  To 
achieve  a  target  accuracy  of  1%  in  this  measure  the  regular  leapfrog  scheme  requires  typically 
15  mesh  points  to  resolve  one  wavelength,  and  about  45  to  meet  the  more  stringent  target 
of  0.1%.  This  criterion  of  N  points  per  wavelength  will  be  used  frequently  in  what  follows 
to  measure  the  resolving  power  of  an  algorithm.  Note  that  N  =  27t/0. 


2.2.2  Stencil  lb 

The  stencil  in  Figure  1(b)  leads  to  a  scheme  that  is  the  simplest  member  of  the  class 
considered  by  Iserles  [1].  On  this  stencil  «|  is  approximated  using  the  average  evaluated 
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over  the  two  time-like  legs,  and  u,  using  the  single  space-like  leg.  After  collecting  terms, 
we  obtain  the  scheme 

= (1  -  2*')  («”+i  -  <)  •  (5) 

The  dispersion  relationship  for  this  scheme  is 

sin  —  912)  -1-  (1  —  2v)  sin  9/2  =  0. 


For  stability,  we  see  that  we  must  have 


-1  <  1  -  2»/  <  1, 


implying  0  <  »/  <  1.  '  This  stencil,  as  we  might  have  guessed,  is  appropriate  only  for 
right-going  waves.  For  left-going  waves  the  mirror  image  of  this  stencil  is  needed. 

The  leading  term  of  the  phase  error  is 

£p^\^v(l-v){l-2v) 

Target  errors  of  1.0%  or  0.1%,  are  met  for  all  i/  by  a  resolution  of  A  >  6  or  >  18 

respectively.  This  scheme  is  roughly  three  times  as  economical  as  regular  leapfrog. 


2.2.3  Stencil  Ic 

Further  improvements  are  possible  by  creating  higher-order  schemes.  On  the  stencil  of 
Figure  1(c)  the  time  derivative  can  be  approximated  as  previously,  but  the  space  derivative 
may  have  weight  k  given  to  the  inner  interval  and  weights  ^(1  -  k)  to  each  outer  interval. 
The  resulting  scheme  is 

“i«  -  «r'  =  (1  +  ■'<>  -  +  ■'(1  -  k)  («;«  (6) 

with  dispersion  relationship 

sin  {<f>  —  9/2)  -1-  (1  -I-  v{\  -  3A:))sin0/2  —  v{\  -  l:)sin3^/2  =  0 


This  scheme  has  a  leading  phase  error  of 
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If  this  is  eliminated  by  choosing 


k  = 


2  3i/  5  31:1 

~  T~  ~2\' 

5  -i-  Oi/  —  Au^ 


we  obtain  the  scheme 

=  -i(l+<')(2-./)(l-2«')  («?+,  -  «")  +  iKl-‘')(l-2i')  («?+2  -  <-.)  ,  (7) 
which  was  also  noted  by  Iserles.  The  dispersion  relationship  can  be  arranged  to  read 
1  -  9/2)  =  {2v  -  l)sin^/2  j^l  +  |t'(l  -  i/)sin^^/2. 


sin  I 


'  Iserles  gives  the  stability  limit  as  0  <  i/  <  j  since  his  stability  definition  encounters  a  technical  objection 
when  V  =  ^.  But  since  the  scheme  reduces  in  that  case  to  the  objection  is  insubstantial. 
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Stability  requires  that  the  RHS  remains  in  the  interval  [-1,1]  and  it  is  clear  that  this  will 
be  so  for  0  <  »/  <  1. 

The  phase  error  of  this  scheme  has  the  leading  term 

ep  =  ^*^(1  -  “  2i/)(2  -  I/) 

The  target  accuracies  of  1.0%  and  0.1%  are  achieved  by  taking  iV  >  4  and  N  >  6  respec¬ 
tively. 


2.2.4  Stencil  Id 


A  different  way  to  achieve  fourth-order  accuracy  is  to  adopt  the  stencil  of  Figure  1(d).  Here 
the  two  alternative  evaluations  of  Ux  must  be  given  equal  weight,  but  the  two  evaluations 
of  ut  from  the  middle  of  the  stencil  can  be  given  weight  ^k,  while  the  evaluations  from  the 
upper  and  lower  legs  each  have  weight  |(1  —  k).  This  leads  to  the  scheme 

(1  -  t)  («J+'  -  «J-»)  =  (1  -  ^  -  2*)  +  (■'-*)  (“"  -  »”«)  •  (8) 

The  dispersion  relationship  can  be  written 

(1  —  k)sm  — ^  +  (2k  +  u  —  l)sin  -  (i/  -  k)sm  =  0. 

The  leading  error  term  is 

26^ 

-—1/(1  -  i/)  [5t/  —  1  —  6fct/] . 

a 

After  removing  this  by  the  choice 

,  5i/- 1 


the  scheme  becomes 


(1  +  i/)  («”+*  -  u”-"*)  =  2il+u)(l-3u)  (uj+i  -  u]^-»)-Kl-3i/)(l-2t/)  [u]  -  .  (9) 

and  the  dispersion  relationship  becomes 

(1  -I-  i/)sin  ^  -  2(1  -f  i/)(l  —  3i/)sin  ^  ^  ~  (1  “  3t/)(l  -  2j/)sin  .  (10) 

The  leading  error  term  is  then 

In  this  case,  the  target  of  1.0%  error  is  met  even  by  =  2,  and  to  achieve  0.1%  requires 
merely  A  >  3. 

To  investigate  stability,  the  dispersion  relation  (lO)can  be  rearranged  as 


tan^/2  =  tan<^/2 


3i/  —  (1  —  2i/)  tan^  <^/2 
3i/*  —  (1  -I- 1/  —  3t/2)  tan^  ^/2 


For  any  real  value  of  tan  0/2,  stability  requires  three  real  values  of  tan<^/2  (otherwise,  as 
discussed  above,  one  of  the  two  complex  roots  is  unstable).  By  considering  the  graph  of 
the  expression  on  the  RHS,  it  is  easy  to  show  that  three  real  roots  are  found  if  and  only 
if  the  singularities  in  the  function  interlace  with  its  zeros.  This  implies  the  stabilty  bound 

0<»/<  i. 
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2.3  Phase  Error  Diagrams 

For  the  four  schemes  discussed  above,  Figure  2.3  shows  the  phase  error  as  defined  in  (4). 
The  improvement  at  each  stage  is  striking. 

2.4  Interpolation  Properties  of  Stencils 

For  linear  advection,  all  of  the  above  schemes  could  b.?  derived  by  using  the  fact  that  the 
solution  is  constant  along  characteristics  dxfdt  =  a.  Referring  to  Figure  3,  the  black  circles 
denote  places  where  the  data  is  available,  and  the  open  circles  places  where  it  can  found 
by  extrapolating  along  a  characteristic.  A  polynomial  can  be  fitted  to  this  combination  of 
values,  and  an  interpolated  value  can  be  read  off  at  the  open  square.  This  value  can  then 
be  extrapolated  to  the  black  square,  where  the  solution  is  required.  The  whole  process 
can  be  thought  of  as  maximum-order  polynomial  interpolation  in  the  coordinate  x  -  at.  It 
is  illuminating  to  compare  the  stencils  on  the  left  with  those  on  the  right.  In  the  latter, 
the  data  within  which  interpolation  must  be  carried  out  are  much  more  compactly  located. 
This  is  the  basic  reason  for  the  improvement  in  going  from  (a)  to  (b)  and  from  (c)  to  (d). 

For  problems  without  constant  coefficients,  the  interpretation  will  not  be  so  simple,  but 
the  advantages  of  the  compact  stencil  will  remain.  A  desire  to  group  the  stencil  tightly 
around  the  characteristic  wiU  lead  to  upwind  biassing. 

2.5  Group  Velocity  and  Boundary  Procedures. 

When  waves  propagate  without  dissipation,  the  concept  of  group  velocity  applies.  This  can 
be  variously  interpreted  as  the  speed  of  a  wave  packet,  or  the  speed  with  which  energy 
propagates.  Linear  first-order  PDFs,  without  source  terms,  are  not  dispersive,  and  so  their 
group  velocity  is  the  same  as  their  phase  speed.  Trefethen  [2]  has  drawn  attention  to  the 
importance  of  group  velocity  for  discretized  PDFs,  which  are  dispersive.  In  particTiIar,  if  the 
group  velocity  should  have  the  opposite  sign  from  the  phase  velocity  at  any  wavenumber, 
such  waves  are  prone  to  be  introduced  into  the  solution  as  spurious  oscillations  by  numerical 
boundary  conditions. 

Iserles  [1]  gave  two  interesting  theorems  relating  to  generadised  leapfrog  schemes  with 
three  time  levels.  The  first  of  these  states  that  if  the  scheme  is  to  be  stable  for  small 
Courant  numbers,  the  gap  between  the  two  time-like  legs  must  comprise  either  zero,  one, 
or  two  space  intervals.  He  then  proved  that  with  zero  or  two  intervals,  the  group  velocity 
would  always  be  in  the  wrong  direction  for  the  shortest  waves  with  6  —  ir.  It  has  the  correct 
sign  for  stencils  with  one  interval. 

The  group  velocity  is  defined  by 

(11) 

Computing  this  expression  for  =  tt  in  each  of  the  four  cases  yields 


_  d<f>Ax 
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Figure  2;  (a)  Phase  error  as  function  of  resolution  for  the  regular  Leapfrog  scheme.  The 
dashed  horizontal  lines  correspond  to  one  part  per  thousand. 


Figure  2:  (b)  Phase  error  as  function  of  resolution  for  the  Upwind  Leapfrog  scheme. 


Figure  2:  (c)  Phase  error  as  function  of  resolution  for  the  Upwind  Leapfrog  scheme  extended 
in  space. 
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Figure  2:  (d)  Phase  error  as  function  of  resolution  for  the  Upwind  Leapfrog  scheme  extended 
in  time.  Note  that  the  Courant  numbers  have  been  halved. 
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Figure  3:  Interpretation  of  the  schemes  as  interpolation  schemes 
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Figure  4:  Group  velocities  of  the  four  schemes  for  short  waves. 
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Figure  2.5  shows  these  short  wave  group  velocities  as  functions  of  Courant  number. 
Cases  (a),  (b),  and  (c)  exemplify  Iserles’  theorem.  Case  (d),  which  is  not  covered  by  the 
theorem,  does  best  of  aU,  the  group  velocity  being  exact  when  u  =  ^,5.  Probably  the  exact 
value  does  not  matter,  since  the  schemes  should  not  be  used  to  model  waves  that  are  so 
underresolved,  but  all  three  of  the  upwind  leapfrog  schemes  display  behaviour  that  leads  to 
no  great  anxiety  regarding  boundary  problems. 


2.6  Second  (Spurious)  Solutions 

A  well-known  danger  of  using  multi-level  schemes  is  that  the  dispersion  relationship  has 
spurious  branches.  Here,  in  the  case  of  stencils  (a),(b),  or  (c),  if  the  phase  angles  {0, 4>)  yield 
a  solution,  then  so  do  the  angles  (0,<p+  rr).  One  of  these  solutions  provides  a  reasonable 
approximation  to  the  continuous  problem;  the  other  does  not.  For  given  0,  the  amplification 
factors  of  the  two  solutions  are  in  the  ratio 

exp  ii<p+  t)  _ 
exp  i{4>) 

so  that  the  second  solution  can  be  expected  to  oscillate  in  time.  Such  a  solution  can  enter 
the  calculation  either  if  (a)  a  poor  starting  procedure  is  used  at  the  first  time-step,  or  (b) 
through  aliasing  caused  by  nonlinearity,  or  (c)  through  the  boundary  procedures.  However, 
in  no  calculation  made  to  date  with  the  upwind  leapfrog  schemes  has  any  contamination 
occured  from  these  sources. 

What  have  been  observed  are  spurious  solutions  due  to  source  terms.  This  danger  can 
easily  be  demonstrated  on  the  model  problem 

Ut  4-  aux  =  su  ( 12) 


whose  exact  solution  is  u{x,t)  =  e*‘u(x  -  at,0).  A  naive  extension  of  (5)  would  be 

“j+i  ~  =  (1  -  2'")  (“j-n  -  “j)  +  -f  Uj+i).  (13) 

The  ansatz  u”  =  g'^expij0  leads  to  a  quadratic  equation  for  the  amplification  factor  g, 
with  roots  51,2  whose  product  is  given  by 


9i92  = 


q 


n— I 


C?,‘ 


=  1 


where  C"  is  the  coefficient  of  u”  in  the  difference  scheme.  Let  gi  be  the  amplification  factor 
of  the  true  solution,  and  g2  that  of  the  spurious  solution.  Suppose  that  s  <  0,  so  that  the 
true  solution  decays.  Then  <  1,  but  </2  >  1,  and  the  spurious  solution  will  grow.  Even 
though  the  mechanism  for  introducing  it  may  be  very  weak,  it  will  eventually  dominate  the 
calculation. 

A  way  to  remove  the  spurious  solution  from  the  model  problem  is  to  eliminate  the  source 
term  by  making  the  substitution 


u  =  e’^v  =>  vt  +  avx  =  0 


(Note  that  this  still  works  if  s  depends  on  x,  and  even  if  u  is  a  vector  and  s  a  matrix.)  Now 
solve  for  v  using  (5)  and  then  translate  back  to  «.  The  result  is 

=  (1  -  21/)  -  u])  .  (14) 


By  construction,  both  roots  now  have  the  same  amplification  factor  (which  is  exact),  and 
the  spurious  root  can  never  come  to  dominate.  Equation  (14)  can  be  rearranged  to  show 
the  source  term  more  explicitly  as 


-  u”-’)  +(1  -e-’‘^‘)«7+,  -  (1  - 


(15) 
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In  this  form,  precise  and  elaborate  instructions  appear  for  evaluating  both  the  time  deriva¬ 
tive  and  the  source  term.  However,  the  far  simpler  scheme  in  which  is  replaced  by 

1  ±  sAt  also  has  the  property  that  =  |j2|,  in  fact 


\9i\  =  Iff2|  = 


'1  -I-  sAf 
1  —  sAi ' 


which  approximates  to  second  order.  This  provides  a  practical  algorithm  ia  which  the 
spurious  solution  is  always  kept  under  control,  provided  sAt  is  not  large. 


2.7  General  Remarks 

The  Upwind  Leapfrog  methods  have  a  number  of  attractive  features  as  starting  points  for  the 
construction  of  practical  wave  propagation  algorithms.  Chief  among  these  is  the  fact  that 
clustering  th?  stencil  around  the  characteristic  enables  high  accuracy  to  be  achieved  with 
a  low  operation  count  in  a  fully  discrete  way.  For  application  to  a  system  of  equations  the 
left-  and  right-going  characteristic  e^-ations  must  be  discretised  on  different  stencils.  The 
overhead  for  doing  this  is  the  cost  of  one  simple  Riemann  solution  per  pair  of  adjacent  points. 
The  second-order  Upwind  Leapfrog  method  has  been  applied  in  this  way  to  the  linear,  but 
non-constant-coefficient  problem  that  arises  from  placing  a  noise  source  in  a  quasi-one- 
dimensional  nozzle  [3].  The  version  described  there  was  later  found  to  suffer  long-time 
instability  due  to  the  second  solution  discussed  above,  but  the  modification  described  has 
allowed  very  long-time  calculations  to  be  made  (Choolwan  Kim,  private  communication). 

Upwinding  normally  pays  dividendr  at  boundaries.  Especially  as  regards  Schemes  (b) 
and  (d),  the  use  of  a  two-point  stencil  ia  space  means  that  no  special  procedures  are  called 
for  at  inflow  or  outflow  boundaries.  Scheme  (c)  should  be  simpler  in  this  regard  than 
other  fourth-order  methods.  Another,  related,  benefit  is  the  ability  to  handle  locally  refined 
meshes  without  generating  reflected  waves. 

Multilevel  methods  require  special  starting  procedures  for  the  first  few  time  steps.  For 
initial- value  problems  a  conventional  high-order  scheme  can  be  used,  but  if  the  solution  is 
excited  from  the  boundaries,  as  with  many  wave  propagation  problems,  no  problem  arises. 

All  of  the  methods  are  conservative,  in  the  sense  that  for  a  scheme  with  /  levels 

-I-  boundary  terms. 

j  j 

because  all  terms  at  intermediate  time  levels  ‘  telescope’.  If  a  conservative  starting  scheme 
is  employed,  or  if  all  perturbations  originate  at  the  boundary,  the  schemes  are  conservative 
in  the  conventional  sense.  However,  without  a  limiting  procedure,  they  cannot  be  used  for 
nonlinear  problems  where  shocks  may  form. 

Scheme  (d)  is  not  suitable  for  steady-state  problems,  because  its  stencil  is  not  broad 
enough  to  give  better  than  second-order  accuracy  in  space.  However,  for  wave  propagation 
it  keeps  its  formal  accuracy  even  when  u  is  vanishingly  small. 

It  is  worth  noting  that  for  stencils  (b)  and  (d),  the  order  of  accuracy  is  unaffected  by  any 
mesh  nonuniformity,  since  onlj  one  mesh  interval  is  involved.  Of  course  regular  timesteps 
must  be  taken. 

These  considerations  provide  strong  motivation  for  attempting  the  extension  to  higher 
dimensions,  which  will  be  the  topic  of  the  remaining  sections. 
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Figure  5:  Possible  stencil  for  linear  advection  in  two  dimensions. 


3  Two-dimensional  Problems 


The  natural  way  to  extend  upwind  leapfrog  schemes  to  advection  in  two  dimensions  turns 
out  not  to  be  very  fruitful.  Consider  the  stencil  in  Figure  5.  This  appears  to  offer  a  suitable 
basis  for  solving 

Ut  +  +  6tij,  =  0 

if  a,  6  >  0.  There  is  a  unique  reversible  discretization  on  this  stencil  that  is  indeed  second- 
order  accurate,  and  also  stable  provided  aAt/Ax,  bAt/Ay  €  [0, 1].  However,  the  scheme  is 
dificult  to  apply  in  situations  in  which  a,  6  are  not  constant.  For  example,  suppose  that  b  is 
smaU  and  frequently  changes  sign,  thereby  switching  the  stencil  between  two  alternatives; 
if  there  is  a  significant  curvature  in  the  y-direction  this  will  make  for  erratic  results. 

Paradoxically,  though,  it  has  proved  possible  to  adapt  the  upwind  leapfrog  schemes 
for  sets  of  coupled  wave  equations  in  more  than  one  dimension.  A  stencil  like  that  in 
Figure  5  is  used  to  discretise  a  bicharacteristic  form  of  the  equations.  Because  the  chosen 
bicharacteristics  are  related  to  the  grid  rather  than  the  solution,  no  switching  occurs. 

The  presentation  will  concentrate  on  generalisations  of  the  scheme  defined  by  stencil  lb. 
Generalisation  of  the  regular  leapfrog  scheme  (stencil  la)  of  course  presents  no  problems. 
A  generaUsation  of  the  fourth-order  scheme  defined  by  stencil  Ic  is  possible,  and  forms  par^ 
of  the  thesis  work  of  my  student  J.  P.  Thomas.  At  present  it  appears  difficult  to  generalise 
the  scheme  defined  by  stencil  Id. 

Before  presenting  the  actual  discretisation,  however,  it  is  neccessary  to  make  some  re¬ 
marks  on  the  general  subject  of  bicharacteristic  equations  and  their  discrete  versions. 


3.1  Bicharacteristic  Equations 


Consider  a  hyperbolic  set  of  partial  differential  equations  in  D  space  dimensions  plus  time 
(ari ,  xj,  •  •  • ,  ar/j;  t),  written  as 


dn  n 


OXd 


(16) 


Premultiply  this  equation  by  h ,  a  left  eigenvector  of  A\  having  eigenvalue  Ai .  This  involves 
no  loss  of  generality,  since  we  can  always  rotate  the  coordinates  to  align  xi  with  any  required 


direction.  The  result  is 


(17) 


This  is  called  a  characteristic  equation.  The  number  of  independent  variables  has  been 
reduced  by  one  compared  with  the  original  form.  The  characteristic  equation  holds  in 
a  D- dimensional  subset,  or  manifold,  of  the  D  +  1-dimensional  space-time.  Derivatives 
normal  to  this  manifold  do  not  appear  in  this  equation.  They  are  allowed  to  be  discontinuous 
across  the  manifold,  which  is  thereby  identified  as  a  wavefront.  Reduction  of  the  dimension 
by  more  than  one  is  not  possible  unless  an  1  can  be  found  that  is  a  left  eigenvector  for  more 
than  one  A. 

Special  significance  is  often  given  to  writing  the  characteristic  equation  in  a  particular 
form.  Derivatives  within  the  characteristic  manifold  can  be  written  along  any  convenient  set 
of  geometric  basis  vectors.  One  element  of  the  basis  set  can  be  chosen  as  the  intersection  of 
the  characteristic  manifold  with  a  neighbouring  characteristic  manifold.  The  characteristic 
equation  then  has  the  appearance  of  an  ordinary  differential  equation  along  that  intersection, 
together  with  “source  terms”,  and  is  called  a  bicharacteristic  equation. 

Numerical  methods  based  on  characteristic  or  bicharacteristic  equations  have  often  been 
proposed  in  the  literature.  Good  surveys  of  early  work  can  be  found  in  the  survey  paper 
of  Chuskin  [3]  and  the  thesis  of  Camarero[4].  More  recent  work  has  been  done  by  [5,6]. 
However,  the  clear  connection  between  physical  interpretation,  differencing  techniques  and 
solution  quality  that  is  found  for  one-dimensional  problems  has  remained  largely  elusive. 
The  methods  work,  but  show  no  striking  advantage.  Reasons  for  this  will  be  discussed  later, 
but  we  next  give  two  examples  of  problems  formulated  as  (bi)characteristic  equations.  Only 
the  two-dimensional  case  is  presented,  but  the  extension  to  three  is  natural. 


3.2  Acoustic  Waves 


The  governing  equations  for  acoustic  waves  in  a  uniform  medium  can  be  written  in  the  form 
(2)  with  unknowns  u  =  (u,v,p)^  corresponding  to  two  velocity  components  and  pressure. 
The  matrices  are 


(18) 


where  p  is  the  density  of  the  medium,  and  a  the  sound  speed.  The  eigenvalues  of  Ai  are 
A  =  ±0,0.  The  eigenvector  corresponding  to  A  =  0  does  not  yield  an  interesting  result,  but 
the  other  left  eigenvectors  are  (po,  0,  ±1)  giving  the  characteristic  equations 


0 

0 

l//> 

0 

0 

0 

II 

0 

0 

0 

,  A2  = 

0 

0 

1/^ 

pa} 

0 

0 

0 

pa} 

0 

(du  du\ 

2dv 

(19) 

=  -/>oV 

dy 

fdu  du\ 

odv 

(20) 

These  are  the  characteristic  equations  for  plane  waves  travelling  in  the  ±x-directions. 
In  this  case,  they  are  also  the  bicharacteristic  equations.  The  bicharacteristic  equations  for 
waves  travelling  in  the  ±y-directions  are 


dp  dp  (dv  dv\ 

+  =  -^Tx 


(21) 
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(22) 


dp  dp  (dv  dv\  2^^ 

For  an  arbitrary  direction  of  travel  the  result  is  found  through  multiplication  by  (pa  cos  ff,  pa  sin  1 ). 


3.3  Elastodynamic  Waves 

Consider  an  isotropic  elastic  medium  undergoing  small  transient  displacements.  To  apply 
the  characteristic  analysis  the  governing  equations  must  be  written  as  a  first-order  system. 
There  are  many  ways  to  do  this,  according  to  the  choice  of  dependent  variables.  Here 
we  choose  u  =  (u,v,p,q,r)  where  again  u,v  are  velocities,  p,q  are  the  normal  stresses  on 
surfaces  of  constant  *,  y,  and  r  is  the  shear  stress.  For  this  problem  we  have 


Ai  =  - 


0 

0 

1/p 

0 

0 

0 

0 

0 

0 

1/p 

0 

0 

0 

0 

1/p 

0 

0 

0 

1/p 

0 

pcj 

0 

0 

0 

0 

,  A2  =  - 

0 

apcl 

0 

0 

0 

apcl 

0 

0 

0 

0 

0 

pci 

0 

0 

0 

0 

pci 

0 

0 

0 

pci 

0 

0 

0 

0 

(23) 


Here,  Ci  is  the  speed  of  propagation  for  longitudinal  waves,  which  bring  about  irrota- 
tional  displacements,  and  Cj  is  the  propagation  speed  for  transverse  waves,  which  produce 
divergence-free  dispacements.  The  parameter  a  is  defined  by  a  =  1  —  2(c2/ci)*,  and  is 
about  1  /3  for  most  metals. 

The  eigenvalues  of  Ai  are  (±ci,±C2,0).  Again  the  zero  value  does  not  give  an  inter¬ 
esting  result;  the  remaining  left  eigenvectors  are  (pcj, 0,  ±1,0,0)  and  (0,pc2,0,0,  ±1).  The 
bicharacteristic  equations  for  waves  propagating  in  the  ±x-diiection6  are 


dp 

dp 

1 

du\ 

2dv 

dr 

(24) 

-  PCX  ^ 

,dt 

dp 

dp 

/ 

'du 

du\ 

2dv 

dr 

(25) 

dt 

f  pci  ( 

^dt 

=  ^^^Yy 

longitudinal  waves,which  are  also  caUed  primary  waves  or  P- 

waves,  and 

dr 

dr 

{  dv 

dv\ 

2du 

dq 

(26) 

di 

-  pC2 

^^^Yx) 

=  ^^Yy 

-  C25- 
dy 

dr 

dr 

/  dv 

dv\ 

2^« 

(27) 

dt 

-^^Yx 

+  PC2 

-^^Yx) 

= 

¥C2-^ 

dy 

(28) 

which  relate  to  the  transverse  waves  (aka  secondary  or  S-waves).  It  is  easy  to  write  out  the 
four  bicharacteristic  equations  for  propagation  in  the  ±y-directions. 

3.4  Discussion  of  Bicharacteristic  Discretisations 

The  numerical  use  of  multidimensional  characteristic  equations  has  been  attended  by  contin¬ 
uing  disagreement  over  such  matters  as  which  (and  even  how  many)  characteristic  equations 
to  use,  and  how  to  discretize  them. 

As  in  the  examples  above,  any  characteristic  equation  is  a  relationship  between  interior 
derivatives  in  some  plane.  Ideally  the  stencil  on  which  it  is  discretized  would  comprise  only 
points  in  that  plane.  Failing  that,  the  points  should  lie  as  close  to  the  plane  as  possible; 
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Figure  6:  Characteristic  planes  for  waves  moving  in  the  positive  i-direction,  and  the  char¬ 
acteristic  cones.  Acoustic  problem  at  left  and  elastodynamic  problem  at  right. 

this  lesson  can  be  learned  from  the  one-dimensional  case.  Bearing  this  in  mind,  we  consider 
the  issues  raised  above. 

It  is  sometimes  pointed  out  that  although  there  are  infinitely  many  bicharacteristic 
equations,  corresponding  to  all  possible  propagation  directions,  only  a  few  of  these  carry 
independent  information.  Indeed,  since  the  bicharacteristic  equations  are  merely  linear 
combinations  of  the  given  PDFs,  it  is  clear  that  not  more  than  three  of  them  can  be 
independent  for  the  acoustic  problem,  or  five  for  the  elastodynamic  problem.  However, 
this  argument  applies  to  the  PDFs  themselves,  acting  at  a  point.  The  discretized  versions 
should  be  related  to  finite  regions  of  space.  It  makes  no  algorithmic  sense  to  ‘add’  the 
pair  of  equations  (19)  and  (20)  and  then  discard  the  terms  that  apparently  ‘cancel’  because 
those  terms  should  be  evaluated,  and  in  an  upwind  discretisation  are  evaluated,  at  different 
places.  Their  difference  is  what  gives  the  bicharacteristic  equations  the  power  to  distinguish 
different  kinds  of  information.  On  a  two-dimensional  Cartesian  grid  there  are,  for  the 
acoustic  and  elastodynamic  equations  respectively,  four  and  eight  bicharacteristic  equations 
that  can  be  naturally  discretised  on  compact  stencils.  All  of  these  carry,  in  the  above  sense, 
independent  information,  and  so  ought  only  to  be  discretised  using  points  that  cluster  round 
the  characteristic  plane  concerned. 

The  problem  then  is  that  there  are  more  bicharacteristic  equations  than  unknowns;  four 
versus  three  in  the  acoustic  problem,  and  eight  versus  five  for  elastodynamics.  This  might 
be  called  the  counting  problem.  In  the  existing  literature  it  is  resolved  in  one  of  two  ways. 
One  is  to  consider  a  reduced  number  of  bicharacteristic  equations.  For  example,  in  the 
acoustic  problem  three  bicharacteristics  might  be  considered,  spaced  at  120°.  However,  the 
derivatives  in  these  surfaces  then  have  to  be  obtained  by  interpolation,  using  data  that  is  not 
the  most  appropriate.  Alternatively,  ail  of  the  bicharacteristic  equations  may  be  employed, 
but  solved  in  some  least-squares  sense;  again  this  results  in  the  muddling  of  information 
that  ought  to  remain  distinct. 

The  conflict  can  be  resolved,  in  the  present  context,  by  storing  the  unknowns  in  a 
staggered  manner.  Fxactly  how  this  is  to  be  done  depends  on  the  way  that  different 
variables  are  coupled  in  the  governing  equations,  and  will  be  described  for  specific  examples 
in  the  next  section. 
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Figure  7:  The  storage  of  unknowns  on  a  square  grid,  and  the  stencil  used  to  update  a  u,p 
point  in  the  second-order  scheme.  The  shaded  plane  shows  a  wave  moving  in  the  positive 
i-direction. 

4  Staggered  Variables 

4.1  The  Acoustic  Problem 

Consider  the  square  grid,  with  uniform  spacing  h  in  both  directions,  shown  on  the  left  of 
Figure  7.  The  x-  and  y-velocity  components  are  stored  on  the  mesh  lines  x,p=constant, 
as  is  commonly  done  for  incompressible  flows,  to  achieve  compact,  centered  stencils  for 
aU  terms.  Here,  however,  pressure  is  stored  at  both  velocity  locations  rather  than  in  the 
cell  centers,  because  the  terms  that  appear  in  any  one  bicharacteristic  equation  must  be 
colocated  on  an  appropriate  stencil. 

How  this  is  achieved  is  shown  in  the  right  half  of  the  figure.  The  stencil  for  the  equation 
describing  waves  moving  in  the  positive  x-direction  comprises  the  points  0,  NW,  SW,  and 
W  at  time-level  n,  together  with  OF  at  level  n  -f  1  and  WP  at  level  n  —  1.  At  the  points 
OF,  O,  W  and  WP  the  variables  u,p  are  stored.  At  NW  and  SW  we  have  v  and  p.  Unique 
reversible  discretizations  for  (19)  and  (20)  are 

2^iiPoF  -  Po)  +  pai^OF  -  uo)]  =  -  Pwp)  +  pa{uw  -  uwp)] 

[(Po  -pw)-  pa  divow)] 

^[(poF-po)-pa(«oF-«o)]  =  --^[{PE  -  Pep)- pa{uE -UEp)] 

[{PE  -  Po)  +  P«  div^oj 

where  divow  denotes  the  ‘discrete  divergence’  uq  —  «w  +  watw  —  vsw-  This  provides 
two  equations  for  the  two  unknown  values  at  OF,  and  eliminates  the  ‘counting  problem’ 
previously  2dluded  to.  Employing  a  more  concise  notation,  with  denoting  a  forward  time 
difference  and  6~  a  backward  time  difference,  and  defining  a  Courant  number  v  =  {aAt)lh, 
this  yields  update  formulae 


i*po  =  ^ 

-f  t/(pw  —  2po  +  PS  -  po(divow  +  divso)) 


pa 

2  V''  Pw  +  ^~Pe)  +  -  ^~«w) 


(29) 
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(a)  (b) 


Figure  8:  Staggered  storage  for  (a)  electricaUy,  and  (b)  magnetically  polarised  electromag¬ 
netic  waves. 


S'^uo  =  —-{6~uw  +  6 

-i-j/(divovv  -  divEo  +  (pw  -  PE)/{pa))  (30) 

Computing  time  is  saved  by  storing  the  terms  that  appear  in  the  first  line  of  each  expression 
from  the  previous  timestep.  The  operations  that  remain  are  those  involved  in  a  first- 
order  upwind  scheme,  except  that  the  appearance  of  (div)  rather  than  just  dxU  imparts  an 
authentically  multidimensional  flavour. 

The  scheme  is  completed  by  adding  a  similar  update  step  for  the  points  where  v,p  are 
stored. 

4.2  Electromagnetics 

In  two  dimensions,  this  problem  has  no  independent  interest,  since  Maxwell’s  equations  have 
the  same  structure  as  the  acoustic  equations.  All  that  is  needed  is  to  distinguish  between 
the  electrically  and  magnetically  polarised  cases,  for  which  the  respective  storage  strategies 
are  shown  in  Figure  8 

4.3  Elastodynamics 

The  storage  here  follows  from  the  fact  that  the  bicharacteristic  equations  in  the  x-  and  in¬ 
directions  each  contain  on  their  right-hand  sides  four  of  the  five  unknowns,  thereby  dictating 
the  storage  shown  in  Figure  9.  The  discretisation  is  then  straightforward. 

4.4  Boundary  Procedures 

Consider  any  ceU  of  which  one  edge  lies  on  a  boundary.  Typically,  some  even  number  2k  of 
variables  are  stored  on  that  edge.  To  update  them,  there  are  k  bicharacteristic  equations 
that  can  be  written  for  the  waves  that  reach  the  boundary  from  the  cell.  These  have  to 
be  supplemented  by  k  conditions  on  the  boundary  itself.  In  the  acoustic  example  treated 
below,  one  boundary  is  a  wall  moving  in  the  x-direction  in  a  prescribed  manner,  so  that  u 
is  given.  The  combination  p  —  pau  is  found  from  the  incoming  bicharacteristic,  so  that  p  is 
determined.  The  velocity  v  parallel  to  the  waU  is  not  required;  in  methods  with  colocated 
variables  some  extra  condition  must  be  introduced  to  find  it. 


«e)  +  2^(^  PE  -6  pw) 


16 


Figure  9:  Staggered  storage  for  elastodynamic  waves. 

At  open  boundaries  there  are  difficulties  that  go  beyond  numerical  considerations,,  since 
for  general  wave  propagation  problems  no  exact  local  boundary  conditions  exist.  A  crude 
expedient,  adopted  in  the  later  examples,  is  to  set  to  zero,  at  the  boundary,  aO  wavestrengths 
from  outside  the  boundary.  An  improved  procedure  is  under  investigation. 


4.5  Von  Neumann  Analysis 

A  von  Neumann  analysis  of  the  schemes  is  possible,  but  complicated  by  the  double  storage 
of  unknowns.  The  acoustic  case  is  presented  for  illustration,  with  the  values  of  p  that  are 
collocated  with  v  renamed  as  q.  Then  one  can  make  the  trial  solution 


9  j  =^|^pog 

where  j,k  are  counting  indices  in  the  i,y-directions.  Inserting  this  into  (19)  yields 

fD  1  IT\  \  ^  /  •  2^\  •  1/  • 

\P  +  U)  sin— COS  — COS— — (»/  — sin  — )sin—  -V'i'Sin-^  =  0.  (31) 

Similar  relationships  can  be  obtained  from  the  other  three  equations,  and  the  condition 
that  all  four  admit  a  non-trivial  solution  for  P,U,Q,V  leads,  after  a  little  rearrangement, 
to  the  determinant 


sin  i  cos  ^  cos  ^  (i/-sin*  sin  ^ 
(i/— sin*  t)  sin  ^  sin  4  cos  4  cos  ^ 


sin  f  cos  ^  cos  ^  (i/-sin*  ^)sin  ^ 
(v— sin*  ^)sin  ^  sin  ^  cos  ^  cos  ^ 


This  is  not  very  tractable  except  in  the  special  case  u  =  \/2  when  the  dramatic  simpli¬ 
fication 

2  i  1  2 

COS  <p  =  COS  - cos  -r 

2  2 

takes  place.  We  can  use  this  result  to  study  the  propagation  of  a  plane  inclined  wave  with 
9x  =  6  cos  a.  By  =  0sina.  We  obtaun 

<i>  \  .  2n 
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Figure  10:  Pressure  contours  due  to  a  wall-mounted  oscillating  piston. 


showing  that  phase  speeds  better  than  1.0%  or  0.1%  are  achieved  for  all  wave  directions 
with  iV  >  5  or  iV  >  15.  Thus,  the  two-dimensional  method  is  isotropic  to  within  the  phase 
error  of  the  one- dimensional  method. 


5  Numerical  Experiments 


The  first  results  relate  to  the  acoustic  equations  solved  in  axisymmetric  coordinates.  This 
requires  the  addition  of  a  term  v/y  to  the  divergence  terms  in  (29,30).  The  upper,  lower, 
and  right  boundaries  are  open.  On  the  left  boundary  the  fluid  velocity  component  u  is 
prescribed.  For  cell  indices  45  through  56  its  value  at  the  nth  time  level  is  set  to 


u 


n 


=  sin 


2iri/n 
N  ■ 


Elsewhere  it  is  zero.  Thus  the  boundary  conditions  simulate  an  oscillating  piston  mounted 
on  a  plane  wall.  The  parameter  N  corresponds  to  the  number  of  mesh  points  per  wavelength. 
Figure  10  shows  a  snapshot  of  pressure  contours  obtained  on  a  100  x  100  grid  with  iV  =  6 
and  V  =  0.45. 

For  this  particular  problem  there  is  an  analytical  solution  for  the  envelope  of  the  press- 
sure  field  on  the  axis  of  symmetry  [7].  Figure(ll)  compares  this  with  the  numerical  pre¬ 
diction  on  a  larger  grid  (400  x  400)  grid  with  iV  =  8  and  v  =  0.45.  It  can  be  seen  that 
the  waves  travel  over  50  wavelengths  without  any  discernable  loss  of  energy.  However, 
careful  inspection  reveals  that  the  amplitude  of  the  numerical  waves  is  slightly  too  large  in 
certain  regions  (x  ~  0.15,0.45,0.80).  This  is  due  the  numerical  wavespeed  being  slightly 
anisotropic.  Waves  that  reach  the  same  place  by  different  routes  are  therefore  slighly  out  of 
phase,  and  their  mutual  interference  is  not  correct.  Solutions  computed  woth  N  A  (not 
shown)  exhibit  pronounced  spurious  interference  patterns. 
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Figure  11:  Pressure  distribution  on  the  axis  of  symmetry  due  to  a  wall-mounted  oscillating 
piston. 

The  second  test  problem  is  from  elastodynamics.  The  upper  boundary  is  a  free  surface, 
on  which  q  =  t  =  0,  except  at  a  single  point  (upper  left  corner)  where  q  is  specified  as  a 
point  loading; 

9(0  =  -^(^  -  <*)  exp  -2  -  exp  -2  j  . 

This  loading  is  taken  from  [8],  and  simulates  a  positive  impulsive  load,  followed  by  a  negative 
impulsive  load.  The  case  with  truly  impulsive  loading  is  called  Lamb’s  problem.  The  left 
boundary  is  a  plane  of  symmetry,  and  the  right  and  lower  boundaries  are  open.  This  is 
harder  than  the  acoustic  problem  because  of  the  two  dilFerent  wavespeeds.  The  slower  S- 
waves  (transverse  waves)  have  a  shorter  spatial  wavelength  and  so  the  mesh  requirement 
is  driven  by  the  need  to  resolve  them.  The  faster  P-waves  (longitudinal  waves)  then  have 
about  twice  as  many  points  per  wavelength  as  they  would  require  by  themselves. 

Figure  12  shows  the  various  waves  that  are  generated,  through  contours  of  the  quantity 
p  +  q.  The  mesh  size  is  200  x  200,  which  resolves  the  S-wave  over  about  10  points,  and  the 
P-wave  over  about  20.  The  head  wave  is  an  S-wave  generated  by  interaction  of  the  P-wave 
with  the  surface,  and  the  Rayleigh  wave  is  a  feature  that  stays  close  to  the  surface  and 
travels  slightly  slower  than  the  S-wave.  In  this  simulation,  it  has  not  yet  detached  itself 
from  the  S-wave.  More  discussion  of  Lamb’s  problem,  including  an  analytical  solution  of 
the  impulsive  loading  case,  can  be  found  in  texts  on  elastodynamics,  such  as  [9]. 

6  Three-dimensional  Applications 

The  method  applies  naturally  in  three  dimensions;  the  correct  storage  strategies  are  easily 
determined.  In  various  cases,  the  quantities  stored  on  a  face  x  =  const  are;  for  acoustics 
(«,p),  for  electromagnetics  {By,  Bz,  Ey,  Eg),  and  for  elastodynamics  {u,v,w,(rxxt(^xy,<fxz)- 
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Figure  12:  Contours  of  the  sum  of  normal  stresses  for  a  smoothed  version  of  Lamb’s  problem. 

Calculations  of  three-dimensional  electromagnetic  problems  have  been  presented  by  Nguyen 
and  Roe  [10]. 

A  central  leapfrog  scheme  due  to  Yee  [11]  has  been  popular  in  electromagnetics,  which 
also  uses  staggered  storage.  The  strategy  in  that  scheme  is  motivated  by  the  wish  to  achieve 
a  compact  discretisation  for  each  component  of  the  curl  operator,  centered  where  it  must 
be  used  to  create  a  time  derivative.  On  a  cubical  array,  the  center  of  each  face  contains  the 
normal  component  of  the  magnetic  field,  and  the  center  of  each  edge  contains  the  parallel 
component  of  the  electric  field.  Accounting  for  the  way  that  edges  and  faces  are  shared 
reveals  that  each  cube  is  the  site  of  six  pieces  of  information.  The  present  scheme  stores 
twelve.  This  is  outweighed  by  the  ability  of  the  present  scheme  to  take  larger  timesteps  (by 
a  factor  £>2 ,  where  D  is  the  number  of  space  dimensions),  and  to  work  with  coarser  grids 
(by  a  factor  2^).  These  are  the  gains  and  losses  of  the  characteristic  formulation.  However, 
a  complete  comparison  remains  to  be  done. 

7  Conclusions  and  Future  Work 

Building  on  the  one-dimensional  foundations  laid  in  [1],  we  have  presented  a  new  type  of 
algorithm  for  linear  wave  propagation.  The  extension  to  higher  dimensions  requires  the 
resolution  of  a  counting  problem  common  to  many  numerical  exploitations  of  the  bicharac¬ 
teristic  equations.  This  is  done  here  by  adopting  a  particular  form  of  staggered  storage. 

Future  publications  will  document  the  application  of  the  method  to  very  long-time 
computations  in  the  presence  of  source  terms,  and  on  dynamically  adaptive  grids. 
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